Multiple gamma rhythms carry distinct spatial frequency information in primary visual cortex

Gamma rhythms in many brain regions, including the primary visual cortex (V1), are thought to play a role in information processing. Here, we report a surprising finding of 3 narrowband gamma rhythms in V1 that processed distinct spatial frequency (SF) signals and had different neural origins. The low gamma (LG; 25 to 40 Hz) rhythm was generated at the V1 superficial layer and preferred a higher SF compared with spike activity, whereas both the medium gamma (MG; 40 to 65 Hz), generated at the cortical level, and the high gamma HG; (65 to 85 Hz), originated precortically, preferred lower SF information. Furthermore, compared with the rates of spike activity, the powers of the 3 gammas had better performance in discriminating the edge and surface of simple objects. These findings suggest that gamma rhythms reflect the neural dynamics of neural circuitries that process different SF information in the visual system, which may be crucial for multiplexing SF information and synchronizing different features of an object.

In the visual cortex, visual stimuli with a high degree of spatial frequency (SF) information such as bars [18][19][20][21] or contours [14][15][16] have been found to induce synchronization in the gamma band. However, more recent studies have also demonstrated that stimuli with low SF information, such as large uniform surfaces, also induce strong gamma synchronous responses [22][23][24]. The finding that both high and low SFs can generate a gamma rhythm is puzzling because high SF information, such as edges [25,26], is thought to be processed by different a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 neural mechanisms than low SF components, such as surfaces [27][28][29][30][31]. There are 2 possible explanations for this contradiction. The first is that a single gamma rhythm can integrate high and low SF information. The second is that distinct gamma rhythms carry different aspects of SF information. The second hypothesis is plausible because distinct narrowband gamma oscillations have been found in cats, primates, and humans [9,20,[32][33][34][35], among other species. However, the relationship between SF selectivity and specific gamma frequency bands remains poorly understood [36,37].
In this study, we examined how SF information is encoded through gamma band neural activity. Surprisingly, we found 3 distinct narrowband gamma rhythms in the visual cortex with different SF preferences and selectivity. We investigated the response properties and neural origins of these gamma rhythms.

Preparation of anesthetized cats
All procedures were in accordance with the National Institutes of Health Guidelines, and the research protocol was approved by the Biological Research Ethics Committee of Beijing Normal University (ID: IACUC(BNU)-NKLCNL 2017-06) and Fudan University (ID: 2021JS0087). Twelve hours before the experiment, cats were injected with dexamethasone (0.4 mg/kg, subcutaneously). The cats (weighing 2 to 4.5 kg each) were initially anesthetized with isoflurane (5% concentration) and injected with atropine sulfate (0.05 mg/kg, subcutaneously). After we conducted a tracheotomy, the animals were artificially ventilated. During recordings, anesthesia and paralysis were maintained with propofol (2 to 6 mg/kg/h) and vecuronium bromide (0.1 mg/kg/h), respectively. The end tidal CO2 was maintained at 3.5% to 4%, and the body temperature was set at 36.5˚C to 37.5˚C. The eyes of the cats were treated with 1% atropine sulfate solution to dilate the pupils, and they were fitted with appropriate contact lenses with artificial pupils (1-mm diameter) to restrict the amount of luminance and ensure focus on a tangent screen. Ventilation pressure, heart rate, electrocardiogram, and blood oxygen were monitored continuously throughout the experiment. To maintain the health of the cats prior to the experiment, they were injected with dexamethasone antibiotic every day (0.4 mg/ kg) and with atropine sulfate (0.05 mg/kg) every other day. One craniotomy (10 to 15 mm in diameter) was performed over the primary visual cortex (VI; centered at P0-2 and L0, mainly in Brodmann's area 18 (A18)), and another (approximately 6 mm in diameter) was performed over the lateral geniculate nucleus (LGN) (centered at L8, A9).

Electrophysiological recordings
Multielectrode arrays (Utah Array, BAU : PleaseprovidethemanufacturerlocationforBlackrockMicrosys lackrock Microsystems, USA) were used to record signals in A18 from the horizontal direction (parallel to the cortical surface) in 12 cats. The arrays consisted of 6 � 8 or 10 � 10 grids of microelectrodes (1-mm electrode length and 400-micron electrode separation). They were inserted 0.5 to 0.6 mm into the cortex using a pneumatic insertion device. Multielectrode linear arrays (Blackrock Microsystems) were used to record signals in area 18 from the vertical direction (perpendicular to the cortical surface) in another 3 cats. The linear arrays consisted of 24 microelectrodes (spaced 0.1 mm apart, each 0.015 mm in diameter, PAU : PleaseprovidethemanufacturerlocationforPlexoninthesentenceThelineararraysconsi lexon, USA). They were inserted 1 to 1.2 mm into the cortex using an automatic controlled insertion device. The multielectrode linear arrays (Plexon) were also used to record signals in the LGN of 6 cats (in 4 of these, recordings were simultaneously conducted in the V1 using Utah arrays). They were inserted 13 to 14 mm into the cortex using an automatic controlled insertion device. To prevent pulsation, the arrays and the exposed cortex were covered with 1.5% to 2% agar. Signals were amplified using a Cerebus 128-channel system (Blackrock Microsystems). All local field potential (LFP) signals were sampled at 500 Hz per channel with a wide-band front-end filter (0.3 to 500 Hz). LFP data were postprocessed by removing channels that were deemed nonfunctional because of broken electrodes or noise.

Visual stimulation
Visual stimuli were generated using a computer with a Leadtek GeForce 6800 video card and displayed on a cathode-ray tube (CRT) monitor (Sony CPD-G520, mean luminance of 32 cd/ m 2 , 1024 � 768 resolution, refresh rate 100 Hz) placed 40 cm away from the animal's eyes. Luminance nonlinearities were corrected using the luminance meter (photoresearch PR670). Stimuli were presented only to the eye contralateral to the craniotomy. For each experiment, we first mapped the receptive fields (RFs) of the multiunit activity (MUA) responses using short instances of sparse visual noise and then conducted an experiment with drifting gratings and large squares with uniform black/white surface luminance.

RF mapping
After manually mapping the RFs of the recorded channels, we used sparse visual noise to precisely locate the RF center. Specifically, a white or black square (0.5 to 2 degrees of visual angle) was flashed on a gray background at different positions in a pseudorandom sequence (usually on a 13 × 13 or 11 × 11 sample grid). The stimuli were presented with an effective frame rate of 100 Hz so that each sparse noise image appeared for 40 ms. Each square was presented 30 to 60 times. The sequence was cut into small segments based on the trial length. This process enabled us to generate a two-dimensional map for each channel. The averaged responses from the x-axis and y-axis of each map were fitted with a one-dimensional Gaussian function to estimate the center position and radius of the RF (σ of the Gaussian function).

Gratings and uniform stimuli
WAU : PleasecheckwhethertheeditstothesentenceWefirstpresenteddriftinggratingswith:::arecorrect; and e first presented drifting gratings with different SFs (approximately 0.01 to 2.5 cycles per degree, 10 conditions) and directions (−18˚to 180˚, 18 conditions). The stimulus contrast was 90%, and the temporal frequency (TF) was 4 cycles per second. Afterward, all drifting gratings were enlarged (approximately 38˚) to cover the RFs of all recording sites, and the stimuli were shown only to the contralateral eye. The temporal sequence of the visual experiment was prestimulus (0.4 seconds), on stimulus (2 seconds), and off stimulus (0.4 seconds). In 2 cats, a large uniform square stimulus (20˚for one cat and 25˚for the other) (black or white against a gray background) was then presented at different positions, so that the center of the RF for that site fell into the border (edge) or the center (surface) of the stimulus.

Laminar alignment
To align the different probe placements in terms of depth, we used the laminar pattern of the MUA responses to the visual stimuli with different SFs. Because the thickness of the cortex and verticality of the probe differed between the probe placements, we used the relative depth (ReD) [38,39] to represent the cortical position of each recording site. The ReD is the normalized cortical depth, ranging from 0 to 1. Three variables were used to calculate the ReD. First, the recording site with the lowest SF selectivity was considered to be in the middle of the input layer (Cha1) and was assigned a ReD value of 0.5. Second, the cortical surface was considered to be located 0.050 mm above the uppermost channel (Cha2) with significant visually driven spiking responses (signal-to-noise ratio [SNR] > 3 and the 3 consecutive channels below Cha2 were also expected to meet the criterion of SNR > 3) and was assigned a ReD value of 0. Third, the boundary separating the gray and white matter was considered to be located 0.050 mm below the deepest channel (Cha3) with significant visually driven spiking responses (SNR > 3 and the 3 consecutive channels above Cha3 were also expected to meet the criterion of SNR > 3) and was assigned a ReD value of 1. After the 3 signatures (Cha1, Cha2, and Cha3) had been determined, they were used as references when calculating the ReD for each channel. ReD values of 0.375 to 0.625 were assigned as the input layer of A18 in the cat V1 according to the procedure used in previous studies [40,41].

Power spectrum analysis
For each trial, the power spectrum of the LFP response over the period 300 to 2,000 ms after stimulus onset (yielding the stimulus power spectrum) was estimated using the multitaper method [42,43] (time-bandwidth product, 3; tapers, 5; Chronux toolbox (http://chronux.org/ )), which was implemented using custom software written in MATLAB. Essentially, the multitaper method attempts to reduce the variance of spectral estimates by premultiplying the data with several orthogonal tapers known as Slepian functions. The frequency decomposition of multitapered data segments therefore provides a set of independent spectral estimates that, once averaged, provides an ensemble estimate that is more reliable for noisy data.

Coherence analysis
We evaluated spike-field coherence (SFC) by calculating the coherency Cxy between different sites (x and y) as the cross-spectra between signals in x and y (Sxy), normalized by the geometric mean of their autospectra (Sxx and Syy) [8]. This was estimated using the multitaper method [44] (time-bandwidth product, 3; tapers, 5; Chronux toolbox (http://chronux.org/)) and implemented using custom software written in MATLAB. To rule out the possibility that the effect of the SFC was solely the effect of the gamma power, we conducted a control analysis in which we calculated the SFC between spike trains, shuffled relative to the LFP signals for each trial within each stimulus condition.

Model fitting and evaluation
To quantify the responses of the different gamma components (GCs) in terms of different SFs, we used a two-dimensional descriptive model to fit the power spectrum. One dimension was the oscillatory frequency, and the other dimension was the SF.
We used a two-dimensional descriptive model (1 to 5) to fit the power spectrums of the LFP and MUA for different SFs. The LFP power in the gamma range could be modeled by a baseline and the weighted sum of the 3 frequency components (1). The frequency profile of the baseline was modeled as a 1/f function that decreased monotonically according to frequency (2), and the TF profiles of the 3 components were all modeled as Gaussian functions (3). The peak frequency μ of each GC was modeled as sigmoid function (4). The SF profile for the power W of the GC was modeled as the difference between the Gaussian function and the sigmoid function for the peak frequency (5). The power spectrums in the edge or surface condition were fitted using formula (1), and the SF tuning curve for the MUA and baseline power in the 3 gamma frequency bands were fitted by formula (5). The goodness of fit was defined as formula (6). Only the sites with a high goodness of fit (larger than 0.8) were analyzed: We fitted power spectrums with the descriptive model (described by formula 1 with details in formulas 2 to 5) by minimizing the mean square error (MAU : PleasenotethattheabbreviationMSEha SE) between the model prediction and the data (power spectrum) using the fmincon function in MATLAB. The parameters that depend on SF include μ i and W i , and the parameters independent of SF include a, b, c, k, σ i , SF 0 , n i , K 0,i , K i , A 1,i , A 2,i , σ 1,i , and σ 2,i (formulas 4 to 5). Instead of determining the number of GCs before model fitting, we fixed the number of GCs at 3 for model fitting to spectrums from all recorded sites. After model fitting, we counted the number of GCs required to explain the power spectrums from each site based on the SNR of each fitted GC. The SNR of each GC was defined as the max weight of the GC (formula 5) divided by the mean power in the corresponding gamma band under the blank condition. If the SNR of a GC was larger than 3, we counted this as a required GC, and otherwise it was considered a negligible GC. In this way, we were able to determine whether a power spectrum had 0, 1, 2, or 3 GCs.

Decoding the edge and surface of a square
We built a logistic regression model to decode the location within a stimulus (edge or surface) that drove a single site (or group of sites) in a given trial, either by gamma power or mean firing rate. The population response matrix was defined as X, X = (X 1 , . . ., X 3N ) T , where X i was a vector (m by 1, m was the total number of trials) of the neural data (either gamma power or MUA response), and the gamma power was the power of the j th recording site in the low gamma (LG) (i = 3j), medium gamma (MG) (i = 3j + 1), and high gamma (HG) (i = 3j + 2) rhythms, respectively. N was the total number of sites included for decoding. As a control, we built a logistic decoder based on the MUA responses using the same procedure. To keep the size of the response matrix (X) constant, we divided the MUA responses from the j th recording site into 3 time bins (if i = 3j, mean rate 300 to 900 ms after the stimulus onset; if i = 3j + 1, mean rate 900 to 1,500 ms after the stimulus onset; and if i = 3j + 2, mean rate 1,500 to 2,000 ms after the stimulus onset). We aggregated all data together for the decoding operation. The units in the X were either all edge response or all surface response. The logistic regression model used to fit the population response to the square location was as follows: Here, W is a vector (3N by 1) of weights. When p � 0.5, the population is predicted to be driven by the stimulus edge; when p < 0.5, the population is predicted to be driven by the stimulus surface. We trained the decoder by optimizing the weight matrix (W) to minimize the loss function (formula 8).
Here, Y is an M by 1 vector representing the location of the stimulus, where 1 is the edge and 0 is the surface. We used an L1 regularization term to minimize the obtained weights and prevent overfitting. We random separate the dataset into a training set (60%), validation set (15%), and test set (25%). For a regularization parameter (λ), we chose the option that gave the lowest MSE in the cross-validation test among 20 lambda values that were equally spaced (logarithmically) between 0.001 and 40. We used the fminunc function in MATLAB to train the decoder. The results shown here ( Fig 3C, S7 Fig) were obtained with the regularization parameter set at 0.1 for all sites.

Site selection
The inclusion criterion for a site was as follows. The firing rate or gamma power of the site at its preferred SF was significantly higher than that of the baseline (i.e., SNR ratio larger than 3). For the MUA responses, the SNR was defined as the standard deviation of the firing rate for visual stimuli divided by that in the blank condition. For LFPs, the SNR was defined as the standard deviation of the gamma power for visual stimuli divided by the standard deviation of the power in the corresponding gamma frequency band in the blank condition.

Statistical analysis
We used an independent sample t test with the Bonferroni correction to compare the SF selectivity across the 3 GCs, gamma power in the edge and surface conditions, the MUA responses in the edge and surface conditions, and the SNRs of the 3 GCs.
We used the bootstrapping procedure to test the influence of location (in the vertical direction in V1) on the peak power of the 3 gamma rhythms. For the i th bootstrap out of 10,000 iterations (bootstraps, i = 1 . . . 10,000), we used the following procedure. We randomly selected 12 recording sessions in the original dataset (recording sessions N = 12) with the same sessions selected more than once allowed. We then computed the mean locations of the peak power for the 3 gamma rhythms (M L for LG, M M for MG, and M H for HG). After 10,000 repeats, we retrieved 10,000 mean values for the locations with peak power. We tested the differences between M L , M M , and M H in a pair-wise manner according to the null hypothesis that X would be no larger than Y. The p-value was calculated as the proportion of the 10,000 repeats in which X was larger than Y.

Results
We used multichannel arrays (96 or 48 channels, spaced 400 um apart, Utah array, Blackrock Microsystems) to record MUA and LFP signals in V1 while presenting large sinusoidal drifting gratings (>38˚of visual angle, see Materials and methods for more details). Fig 1A shows typical LFP waveforms recorded in response to drifting gratings with different SFs. For most recording sites, drifting gratings were able to induce strong gamma band activities under a range of SFs (Fig 1B).

Three distinct gamma rhythms in the visual cortex
Surprisingly, we found that visual stimuli at certain SFs could simultaneously induce 3 narrowband gamma rhythms in V1 (Fig 1B). This was consistent in different cats (Fig 1C). To dissect and quantify the distinct narrowband gamma rhythms, we built a descriptive model to explain the power spectrum at different TFs induced by different SFs (Fig 1C). We assumed that the LFP power spectrum could be modeled by the summation of a baseline and 3 bell-shaped components peaking at different TF ranges: LG (25 to 45 Hz), MG (45 to 65 Hz), and HG (65 to 100 Hz). The frequency profiles of the 3 oscillatory components (GCs) were all modeled as Gaussian functions, and the frequency profile of the baseline was modeled as a 1/f function that decreased monotonically in frequency. We further assumed that the GCs would change with the SF of the visual stimulus. Accordingly, the amplitudes of the 3 oscillatory components were modeled as the difference of Gaussian functions [45][46][47][48], and their peak frequencies were modeled using the Naka-Rushton function [49][50][51][52] (see Materials and methods for the descriptive model and model fitting).
The power spectrums of the LFPs for individual sites could be well reconstructed by this model (see examples of model fit in Fig 1B and 1C). We found that the LFP could be well   Fig 1D), and among these, 34.5% of the sites (n = 250) had 3 strong (SNR > 3) GCs, 37.9% (n = 275) had 2, 21.4% (n = 155) had 1, and 6.2% (n = 45) had no significant gamma rhythms (Fig 1F, S1 Fig). The 3 gamma rhythms could also be observed in the MUA power spectrum (S2 Fig) and were well fitted by this descriptive model (S3 Fig). After successfully dissecting the 3 GCs, we next demonstrated that the 3 gamma rhythms had distinct response properties in V1.

The 3 gamma rhythms had distinct SF selectivity
We found that the power and peak frequency of the 3 gamma rhythms were tuned to different ranges of SF (see Fig 2A and 2B for sites with all 3 GCs, n = 250). This could be quantified according to the significant differences in the cutoff SF (quantified as the largest SF corresponding to 50% of the peak value) and SF selectivity (quantified as 1 minus the ratio between the response at a low SF and the response at an optimal SF) among the 3 GCs. Specifically, the HG was low pass (cutoff SF: 0.09 ± 0.003; SF selectivity: 0.02 ± 0.004; mean ± standard error), whereas the MG and LG preferred relatively higher SFs (cutoff SF: 0. We next examined how the SF tuning properties of the 3 GCs differed from those of spike activity. We found that both the SF selectivity and cutoff SF for spike activity (MUA) were significantly higher than those of HG (p < 0.001, paired t test) and lower than those of LG

PLOS BIOLOGY
(p < 0.001, paired t test) or MG (p < 0.01, paired t test) activity (Fig 2C-2F). This suggests that SF information carried by the mean firing rate at a given recording site has less diversity and a lower spatial resolution than SF information carried by all 3 GCs together. To test whether the diverse SF information we observed in the 3 GCs was also carried by dynamic spike activity [53], we calculated the phase locking strength between spikes and LFPs driven by different stimulus SFs (S5 Fig). We found that the MUA was phase locked to the LFP in the LG range when driven by a stimulus with high SF and phase locked in the HG range when driven by a stimulus with a low SF. In summary, diverse SF information was not only multiplexed by the 3 gamma bands, but also carried by neural dynamic responses that locked to the phases of the 3 GCs.
We also observed different correlations between spike activity (MUA) and the 3 narrowband GCs in terms of their SF tunings: The cutoff SF of MUA was positively correlated to that of LG (r = 0.64, p < 0.001) and MG (r = 0.57, p < 0.001), but the correlation between the cutoff SF of HG and that of MUA was weak (r = −0.21, r = 0.035) (Fig 2E). A similar result was obtained for the relationship between the SF selectivity of gamma activation and MUA (LG, r = 0.44, p < 0.001; MG, r = 0.29, p = 0.0012; HG, r = 12, p = 0.42) (Fig 2F). These results indicate that the 3 GCs may have different neural originations, which we address in a later section. power spectrums in the "RF on surface" and "RF on edge" conditions for a recording site (the blue curve denotes the edge condition, and the red curve denotes the surface condition). B shows a statistical comparison of the response strength in the "RF on surface" and "RF on edge" conditions for the 3 gamma components (33 sites in 2 animals). C shows that the decoding accuracy for the RF position (surface or edge) increased with the number of sites included (the grey curve represents the outcome when the model was trained with gamma power in the LG, MG, and HG bands, while the black curve shows that for the mean firing rate in 3 time bins). D shows the field-field coherence among recording sites with RFs on the square edge (purple) or square surface (blue). The panels A, C, and D are from one cat. The dataset underlying this figure can be found in https://github.com/CleveHan/Gamma-Oscillation. HG, high gamma; LG, low gamma; MG, medium gamma; RF, receptive field. https://doi.org/10.1371/journal.pbio.3001466.g003

Synchronization of the gamma rhythms differentiated SF-related features of a simple object
The results given in the previous sections indicate that 3 distinct GCs carry different SF information about grating patches. We further tested whether different SF information in a static object could also induce the 3 GCs. We presented a black or white square (25 degrees of visual angle, against a gray background) with uniform luminance for 2 seconds in different positions on the screen. The square was presented such that the RF of a site was either on the square's edge, which had a higher SF information, or on its uniform surface (Fig 3A, left panel shows the black square condition), which contained low SF information. We found that LG emerged when the site's RF was on the edge, while MG and HG emerged when the site's RF was on the surface (Fig 3A, middle and right panels). Across the 3 gamma rhythms, the power of LG was significantly higher (p < 0.0001) at the edge ("RF on edge") than at the surface ("RF on surface") of the stimulus, while the power of MG (p < 0.0001) and HG (p < 0.01) was significantly higher for the stimulus surface (Fig 3B). Similar results were found in the white square condition (S6 Fig), but with weaker LFP gamma power (p < 0.0001) compared with the black square condition [22,23,54].
The results in the previous section indicate that the SF resolution of the 3 GCs in a recorded site might be better than that of the mean spike activity. We further tested this idea by decoding the RF positions (on the surface or edge of a simple square) of a site based on the gamma power or spiking activities of that site (see Materials and methods for details). We found that the 3 GCs together could not only be used to predict whether the RF of a site would be on the edge or surface of a simple square, (decoding accuracy was 81.1% ± 0.4%) but also had better decoding accuracy than that using the MUA (decoding accuracy was 66.3% ± 0.8%) from the same recorded site (see Fig 3C for the result from one animal and S7 Fig for that from another  animal). Interestingly, the decoding performance of the 3 GCs from one recording site was even better than the decoding performance based on spike activity from multiple recording sites (Fig 3C) in our dataset. This suggests that the 3 GCs recorded from a single site might be a synchronized signal that pools SF information from spike activity occurring in many sites.
We further calculated the LFP coherence between 2 sites with RFs that were either both on the edge or both on the surface of a square with uniform luminance. We found a coherence peak in the LG band when the RFs of the sites were both on the square edge. By contrast, we found a coherence peak in the HG band when the RFs of the sites were both on the square surface ( Fig 3D, S8 Fig). Interestingly, coherence peaks in the MG band could be seen for both edge and surface conditions. In brief, we have shown that 3 distinct gamma rhythms not only carry and synchronize different SF-related features for simple objects, but that they also discriminate SF-related features better than spiking activity, given a limited number of recording sites.

Cortical distribution and possible origins of the 3 gamma rhythms
We next asked whether the 3 distinct gamma rhythms had same neural generation mechanisms. The correlation between the SF tuning properties (cutoff SF and SF selectivity) of gamma activation and that of the MUA (Fig 2E and 2F) suggest that the 3 GCs may have different neural originations. In this section, we directly measured gamma rhythms at different cortical depths and in the LGN.
We used linear arrays (U-Probe) to simultaneously record from multiple sites across different layers (vertical direction) in V1 (Fig 4A). We found that the power of the 3 gamma rhythms peaked at different cortical depths (see the right side of Fig 4A for an example). There was a tendency for the HG power to peak in the input layer, the MG power to peak near the upper bound of the input layer, and the LG power to peak above the upper bound of the input layer (right side of Fig 4A). The cortical depths of the 3 gamma peaks were significantly different for the mean values from 3 animals (Fig 4B). The average peak of the LG was located in the superficial layer, that of the HG was in the input layer, and that for the MG fell between the other 2 components (p < 0.005 for all pair-wise comparisons) (Fig 4C).
The positions of the power peaks of the HG and MG components in the input layer suggest that they may have subcortical origins. We next asked whether the HG and MG components originate in the visual thalamus. We used different arrays to simultaneously record from V1

PLOS BIOLOGY
(Utah array) and the LGN (U-Probe) (Fig 5A). Fig 5B shows the results of 2 typical simultaneous recordings: 3 gamma rhythms were found in V1 LFP signals (Fig 5B, top row), whereas only HG was found in the LGN MUA signals (Fig 5B, middle row). The SFC between the LGN and V1 (see Materials and methods) further indicated that only HG was generated in the LGN (Fig 5B, bottom row). Among the 84 recorded sites in the LGN, 83.3% of the MUA spectrum could be well fitted by the model (goodness of fit > 0.8, Fig 5C), and all of the well-fitted examples had only 1 GC in the high frequency range (Fig 5D, S9 Fig). The frequency ranges of the HG components were comparable between the LGN and V1 signals in the input layer ( Fig  5E), indicating that the HG signal in V1 is likely inherited from the LGN. We observed a trend in which SF selectivity for HG in V1 was larger than that in the LGN (S10 Fig), as well as a trend in which the cutoff SF for HG in V1 was lower than that in the LGN. These results suggest that cortical inhibition might occur to modify the SF tuning transform from the LGN to V1 [45,46].

Discussion
To our knowledge, this study is the first to report the presence of 3 narrowband gamma rhythms (distinct from the broadband gamma in previous studies [32,36,37,[55][56][57]) with distinct response properties in the V1. The 3 gamma rhythms were distinguished by their oscillatory frequencies: "low gamma" (LG, approximately 25 to 40 Hz), "medium gamma" (MG, approximately 45 to 60 Hz), and "high gamma" (HG, approximately 65 to 90 Hz). The different GCs carried different aspects of visual information in terms of SF. Their distinct neural origins further suggest that gamma rhythms can be independently generated by multiple dynamical systems. Our findings are summarized in Fig 6.

Neural mechanisms for distinct gamma rhythms
An important finding of this study is that the 3 detected gamma rhythms were generated by different neural circuits. Visually induced gamma band activity is generally thought to reflect network properties of excitatory and inhibitory interactions (E-I), which may exist at many stages of visual system (such as in the retina, LGN, and V1). To date, 2 mechanisms have been proposed regarding the generation of gamma rhythms in the visual cortex. Some studies have suggested that gamma is an emergent property of the cortex [58] and that it is likely to be Fig 6. A schematic illustration of our main findings regarding the response properties and generation mechanisms of the 3 gamma rhythms. When the SF information from the natural images (see left) enters the visual system through the retina, to the SF information might be dissected into 3 different frequency channels (right column) in the superficial layer of V1. Signals in 3 colors (red, green, and blue) denote the 3 gamma components (LG, MG, and HG), respectively. HG, high gamma; LG, low gamma; MG, medium gamma; SF, spatial frequency; V1, primary visual cortex.
Our results further show that the MG component was strongest in the V1 input layer although the LG peaked in the V1 output layer (Fig 4). This suggests that MG might be related to local E-I in the input layer [65][66][67][68][69][70][71][72][73], while LG might be related to horizontal connections in the output layer, which has long-distance connections on a larger spatial scale. The distinct mechanisms of LG and MG are consistent with the model prediction from a recent computational study [34]. The involvement of different neural circuits in V1 gamma rhythms also explains the finding that each gamma rhythm has a specific relationship to the spike activity (S5 Fig) and suggests that different gamma bands might be related to different neural functions.

What is the functional role of gamma rhythms?
Our results indicate that both high [14][15][16][18][19][20][21] and low [22][23][24] SFs can induce strong gamma band activity. Specifically, these data represent the first evidence that different GCs encode specific SF information and that they might synchronize different features. More importantly, we found that the SF selectivity of LG was better than that of the mean firing rate. Accordingly, we propose a potential function for GCs in the visual cortex: Multiplexing SFspecific information synchronized [91][92][93][94][95][96][97] in distinct gamma bands. This proposal might not be contradictory to existing theories.
The 2 current hypotheses about gammas are the "binding by synchronization" (BBS) hypothesis [18,[98][99][100][101][102][103][104][105] and the "communication by phase locking" (communication through coherence, CTC) hypothesis [2,106]. The "BBS" hypothesis emphasizes the integration of visual information about an object via gamma synchronization. However, the fact that SF information about a visual stimulus is carried by different GCs suggests gamma rhythms do not integrate visual features together in the visual cortex, but instead integrate them separately in different frequency bands. Through our coherence analysis, we found strong mutual synchronization between LG activities at recording sites with edge inside their RFs and strong mutual synchronization between HG activities at sites with RFs inside the square stimuli (S8 Fig). These results suggest that although SF-specific information is represented in different frequency bands, the gamma rhythms in each frequency band may still contribute to a routingby-synchrony scheme.
The "CTC" hypothesis emphasizes the idea that information might be carried by different phases of a single gamma rhythm and that phase locking might be used for information integration in downstream cortical regions. Our findings suggest that visual information might be carried by different rhythms. Strong coherence between LGN and V1 in the HG band in our study supports the CTC theory. Our study also showed that there are spike-LFP phase locking for the low and medium GCs in V1. According to previous studies on 2 gammas in V1 [33,34], the spatial scale of the communication for LG and MG might be different.
LG generated in superficial layer of V1 might have a relative larger communication scale than that of MG. But the exact forms of communication in these multiple gamma bands in downstream brain regions should be the focus of future theoretical and experimental studies.
There are 2 possible mechanisms for SF dependence in distinct gamma rhythms. The first is that there is generally a negative correlation between brain signals at different TF bands and their preferences for SFs, and the second is that the SF dependences of the 3 gamma rhythms reflect the characteristics of different neural circuits [34,59,66] that generate them. We verified that SF tuning of gamma rhythms is a property of the oscillatory components instead of the TF (S11 and S12 Figs).

How many gamma rhythms are in the visual cortex?
We found 3 distinct narrowband GCs (LG, MG, and HG), although previous studies have only reported 1 or 2 gamma rhythms in the visual cortex of multiple species [9,18,20,23,33,62,79,81,[107][108][109]. This discrepancy is likely due to the fact that gamma rhythms are dependent on the stimulus SF, and that most studies [18,20,79] did not vary the visual stimuli in a range of SFs. Because the 3 GCs had distinct SF tunings (Figs 1 and 2), stimulus conditions that could induce 3 simultaneous and significant GCs were relatively rare (the third row in Fig 1B). Our results also show that 1 or 2 obvious gamma rhythms could be seen under most stimulus conditions, although this could be seen as 3 gamma rhythms if the response patterns of the GCs are considered as a function of SF (Fig 1C). Another possibility why the 3 GCs might not have been reported previously could be the small size of the stimuli used in past studies [82,85]. Larger visual stimuli could induce stronger gamma rhythms compared with smaller visual stimuli [33,34,66]. The use of large visual stimuli in the present study allowed us to locate multiple significant gamma peaks.
Finally, in some high SF conditions, the LG may peak at around 25 Hz. According to the classical division of different frequency bands (beta: 14 to 25 Hz and gamma: 25 to 100 Hz), the rhythm in 25 Hz may be considered a beta rhythm instead of a GC. However, the peak frequency of the LG is a function of SF that is continuously monotonically decreasing (red curve in Fig 2B)  We used the Hilbert transform to filter the data to calculate the phase in each frequency band for different SFs. The firing rate in each LFP phase is shown in the same scale in the top row. The content of the second row is the same as that in the first row but with separate scales. The third row shows the phase at peak firing rate for different frequencies with different SFs. The fourth row shows the phase locking index, defined as (max firing rate − mean firing rate)/ mean firing rate in different frequencies. The bottom row shows the normalized phase locking index (subtraction of the phase locking index in different SFs from the phase locking index in the blank condition). LFP, local field potential; SF, spatial frequency.